%cd 'E:\MATLAB scripts\Blue noise confidence bands\Rice_Realizations\Rice_Realizations'
filelistA = dir('*.txt');
C = [];
SPECTRA = [];
max_time = 20000/1;
for ii = 1:max(size(filelistA));
    if filelistA(ii).isdir ~= true
        fname = filelistA(ii).name;
        fid = fopen(fname);
        out = textscan(fid,'%s %f %f','delimiter',',','headerlines',1);
        fclose(fid);
    end
    datetime = out{1};
    mass = out{2};
    dt1 = datetime(1);%time experiment started
    sec = etime(datevec(datetime),datevec(dt1));
    dt = 1.2;
    maxsec = floor(sec(end));
    t = 1:dt:maxsec;%%%time into run in seconds
    t = t';
    mass_int = interp1(sec,mass,t);
    M = movmedian(mass_int,5);
    Q = (M(2:end) - M(1:end-1))/dt;%Q = mass flux in grams/sec
    Q = [0;Q];
    if ii == 1 | ii == max(size(filelistA));
        disp(mean(Q))
    end
    er = abs((mean(Q) - 0.3577)/0.3577);
    if er <= 0.1;
        t = t(1:max_time);
        Q = Q(1:max_time)*50;
        T = max(t);%Max time
        nt = size(Q,1);%total number of mass measurements
        Qdetrended = Q - mean(Q);%Centering of data
        ne = 1;%number of ergodic realizations from full data set
        nte = floor(nt/ne);
        Qdetrended_e = [];
        for j = 1:ne
            qdetrended_e = Qdetrended(1+j*nte-nte:j*nte);
            Qdetrended_e = [Qdetrended_e qdetrended_e];
        end
        Qdetrended_e = Qdetrended_e';
        %% MTM Spectra
        nw = 4;%number of windows
        nsim = 2;%number of random simulations based on AR(1) model fit, to figure out confidence bands
        t = t(1:nte);
        t = t';
        MTMe = [];
        for j = 1:ne;
            datax = Qdetrended_e(j,:);
            datax = datax';
            [Mspecred,specred,po,fd,fr,tabMC,tabchi,tabtchi,theored,rho]=RedConf(datax,t,dt,nsim,nw);
            MTMe = [MTMe po];
        end
        if ne > 1
            MTMmean = mean(MTMe');%ensemble average data
        end
        if ne == 1
            MTMmean = MTMe';
        end
        period = 1./fd;
        %% Find change pts
        lnperiod = log(period);%convert data over to log space
        lnMTMmean = log(MTMmean);%convert data over to log space

        dp = 0.001;%pick an increment that is fine, relative to data spacing, to interpolate unevenly spaced data onto.
        lnperiod_int = min(lnperiod):dp:max(lnperiod(2:end));%create vector with spacing dp to interpolate onto.
        lnMTMmean_int = interp1(lnperiod(2:end),lnMTMmean(2:end),lnperiod_int);%Do the interpolation
        mean_power = mean(exp(lnMTMmean_int));
        MTMmean = MTMmean./mean_power;

        FC = findchangepts(lnMTMmean_int,'Statistic','linear','MaxNumChanges',2);
        C1 = exp(lnperiod_int(FC(1)));
        C2 = exp(lnperiod_int(FC(2)));
        C = [C;C1];
        SPECTRA = [SPECTRA MTMmean'];

        %% Plot data
        figure(1)
        loglog(period,MTMmean,'k','linewidth',0.5);
        if ii == 1;
            loglog(period,MTMmean,'g','linewidth',0.5);
        end
        if ii == max(size(filelistA));
            loglog(period,MTMmean,'r','linewidth',0.5);
        end
        ylabel('Spectral Power')
        xlabel('period (sec)')
        axis([2.4 6e3 0.003 8])
        hold on
    end 
end 
%cd 'E:\MATLAB scripts\Blue noise confidence bands\Rice_Realizations\Rice_Realizations'
fd = fd(2:end-1);
period = 1./fd;
SPECTRA = SPECTRA(2:end-1,:);
MTMe_sort = sort(SPECTRA,2);
figure(2)
hold on
F = Curl_cb;
F(243,:) = [1 0 1];
%F(129,:) = [0 0 0];
nx = size(SPECTRA,2);
for i = 1:size(MTMe_sort,2)
    s = ceil(i/nx*size(F,1));
    s = F(s,:);
    loglog(period,MTMe_sort(:,i),'color',s,'linewidth',1)
    axis([2 1e4 1e1 1e5])
end
MTMe_mean = median(MTMe_sort,2);
loglog(period,MTMe_mean,'k','linewidth',1);
set(gca, 'YScale', 'log')
set(gca, 'XScale', 'log')
ylabel('Spectral Power')
xlabel('period (sec)')
colormap(F)
b = colorbar;
b.Label.String = 'Fraction of spectra with power less than';
%% Generate 95% Confidence band
B95 = [];
P_int = 0.5:0.05:1;
Percents = 1:1:size(MTMe_sort,2);
Percents = Percents./size(MTMe_sort,2);
for i = 1:max(size(MTMe_sort))
    Power = MTMe_sort(i,:);
    b95 = interp1(Percents,Power,P_int);
    b95 = b95(10);
    B95 = [B95;b95];
end
loglog(period,B95,'m','linewidth',0.25);
axis([2.4 6e3 0.003 8])
disp(size(MTMe_sort,2))
%% Generate data to compare imposed periodicity signals to confidence levels
%cd 'E:\MATLAB scripts\Blue noise confidence bands\Rice_Realizations\Rice_Realizations\Grid'
Gr = load('Grid.dat');
filelistB = dir('*.txt');
Sig_Strength = [];
for jjj = 1:max(size(filelistB));
    if filelistB(jjj).isdir ~= true
        fname = filelistB(jjj).name;
        fid = fopen(fname);
        out = textscan(fid,'%s %f %f','delimiter',',','headerlines',1);
        fclose(fid);
    end
    datetime = out{1};
    mass = out{2};
    dt1 = datetime(1);%time experiment started
    sec = etime(datevec(datetime),datevec(dt1));
    dt = 1.2;
    maxsec = floor(sec(end));
    t = 1:dt:maxsec;%%%time into run in seconds
    t = t';
    mass_int = interp1(sec,mass,t);
    t = t(1:max_time);
    mass_int = mass_int(1:max_time)*50;
    M = movmedian(mass_int,5);
    Q = (M(2:end) - M(1:end-1))/dt;%Q = mass flux in grams/sec
    Q = [0;Q];
    T = max(t);%Max time
    nt = size(M,1);%total number of mass measurements
    Qdetrended = Q - mean(Q);%Centering of data
    ne = 1;%number of ergodic realizations from full data set
    nte = floor(nt/ne);
    Qdetrended_e = [];
    for j = 1:ne
        qdetrended_e = Qdetrended(1+j*nte-nte:j*nte);
        Qdetrended_e = [Qdetrended_e qdetrended_e];
    end
    Qdetrended_e = Qdetrended_e';
    %% MTM Spectra
    nw = 4;%number of windows
    nsim = 4;%number of random simulations based on AR(1) model fit, to figure out confidence bands
    t = t(1:nte);
    t = t';
    MTMe = [];
    for j = 1:ne;
        datax = Qdetrended_e(j,:);
        datax = datax';
        [Mspecred,specred,po,fd,fr,tabMC,tabchi,tabtchi,theored,rho]=RedConf(datax,t,dt,nsim,nw);
        MTMe = [MTMe po];
    end
    if ne > 1
        MTMmean = mean(MTMe');%ensemble average data
    end
    if ne == 1
        MTMmean = MTMe';
    end
    period = 1./fd;
    period = period(2:end-1);
    MTMmean = MTMmean(2:end-1);
    %% logbinned intpolated data
    lnperiod = log(period);%convert data over to log space
    lnMTMmean = log(MTMmean);%convert data over to log space
    dp = 0.001;%pick an increment that is fine, relative to data spacing, to interpolate unevenly spaced data onto.
    lnperiod_int = min(lnperiod):dp:max(lnperiod(2:end));%create vector with spacing dp to interpolate onto.
    lnMTMmean_int = interp1(lnperiod(2:end),lnMTMmean(2:end),lnperiod_int);%Do the interpolation
    mean_power = mean(exp(lnMTMmean_int));
    MTMmean = MTMmean./mean_power;
    n = Gr(jjj,1);
    [val,idx]=min(abs(period-n));
    sig_strength = MTMmean(idx)/B95(idx);
    Sig_Strength = [Sig_Strength;sig_strength];
end
%cd 'E:\MATLAB scripts\Blue noise confidence bands\Rice_Realizations\Rice_Realizations'
figure(3)
x = [0 1];
y = [1 1];
plot(x,y,'k')
hold on
plot(Gr(11:15,2),Sig_Strength(11:15),'s','MarkerFaceColor', [.2288 .1748 0.4758],'MarkerEdgeColor',[0 0 0],'MarkerSize',15)
plot(Gr(6:10,2),Sig_Strength(6:10),'d','MarkerFaceColor', [.6332 .9919 0.2394],'MarkerEdgeColor',[0 0 0],'MarkerSize',15)
plot(Gr(16:20,2),Sig_Strength(16:20),'^','MarkerFaceColor', [.9953 .6534 0.1958],'MarkerEdgeColor',[0 0 0],'MarkerSize',15)
plot(Gr(1:5,2),Sig_Strength(1:5),'o','MarkerFaceColor', [.4796 .0158 0.0106],'MarkerEdgeColor',[0 0 0],'MarkerSize',15)
ylabel('Relative spectral signal amplitude')
xlabel('Input signal amplitude/mean input rate')
%axis([0 1 0 11])
legend('Threshold','Period = 12sec','Period = 100sec','Period = 250sec','Period = 1000sec','Location','northwest')
bb = colorbar;
bb.Label.String = 'Period of input signal';
colormap(turbo)
set(gca,'ColorScale','log')
caxis([10 1000])